function [like,grad]=consumpjointsearchmlogitWCabsorb(bjoint,Y,ut,Xo,grad_4yr,idx1,idx2,sdemog,S,q)
    b  = bjoint(idx1);
    bo = bjoint(idx2);

    flag = Y>0;

    in_white_collar = ismember(Y,[2 4 7 9 12 14 17 19]);
    prev_WC         = ut.sprevs(:,7);

    lambda                    = zeros(size(flag));
    lambda(flag & prev_WC==0) = exp(Xo(flag & prev_WC==0,:)*bo)./(1+exp(Xo(flag & prev_WC==0,:)*bo));
    lambda(flag & prev_WC==1) = 1;

    b2flg    = 2+[1:ut.number2];
    b4sflg   = 2+[1+ut.number2:ut.number2+ut.number4s];
    b4nsflg  = 2+[1+ut.number2+ut.number4s:ut.number2+ut.number4s+ut.number4ns];
    bwptflg  = 2+[1+ut.number2+ut.number4s+ut.number4ns:ut.number2+ut.number4s+ut.number4ns+ut.numberpt];
    bgwptflg = 2+[1+ut.number2+ut.number4s+ut.number4ns+ut.numberpt:ut.number2+ut.number4s+ut.number4ns+ut.numberpt+ut.numbergpt]; 
    bwftflg  = 2+[1+ut.number2+ut.number4s+ut.number4ns+ut.numberpt+ut.numbergpt:ut.number2+ut.number4s+ut.number4ns+ut.numberpt+ut.numbergpt+ut.numberft];
    bgwftflg = 2+[1+ut.number2+ut.number4s+ut.number4ns+ut.numberpt+ut.numbergpt+ut.numberft:ut.number2+ut.number4s+ut.number4ns+ut.numberpt+ut.numbergpt+ut.numberft+ut.numbergft];
    bwcflg   = 2+[1+ut.number2+ut.number4s+ut.number4ns+ut.numberpt+ut.numbergpt+ut.numberft+ut.numbergft:ut.number2+ut.number4s+ut.number4ns+ut.numberpt+ut.numbergpt+ut.numberft+ut.numbergft+ut.numberwc];
    bgwcflg  = 2+[1+ut.number2+ut.number4s+ut.number4ns+ut.numberpt+ut.numbergpt+ut.numberft+ut.numbergft+ut.numberwc:ut.number2+ut.number4s+ut.number4ns+ut.numberpt+ut.numbergpt+ut.numberft+ut.numbergft+ut.numberwc+ut.numbergwc];

    % for debugging
    debug = false;
    if debug==true
        disp('ut.number2');
        disp(ut.number2);
        disp('ut.number4s');
        disp(ut.number4s);
        disp('ut.number4ns');
        disp(ut.number4ns);
        disp('ut.numberpt');
        disp(ut.numberpt);
        disp('ut.numbergpt');
        disp(ut.numbergpt);
        disp('ut.numberft');
        disp(ut.numberft);
        disp('ut.numbergft');
        disp(ut.numbergft);
        disp('ut.numberwc');
        disp(ut.numberwc);
        disp('ut.numbergwc');
        disp(ut.numbergwc);
        disp('size(ut.X2nw)');
        disp(size(ut.X2nw));
        disp('b2flg');
        disp(b2flg);
        disp('b4sflg');
        disp(b4sflg);
        disp('b4nsflg');
        disp(b4nsflg);
        disp('bwptflg');
        disp(bwptflg);
        disp('bwftflg');
        disp(bwftflg);
        disp('bwcflg');
        disp(bwcflg);
        disp('bgwptflg');
        disp(bgwptflg);
        disp('bgwftflg');
        disp(bgwftflg);
        disp('bgwcflg');
        disp(bgwcflg);
        assert(length(b2flg)==32,  'b2 flag wrong');
        assert(length(b4sflg)==32, 'b4s flag wrong');
        assert(length(b4nsflg)==32, 'b4h flag wrong');
        assert(length(bwptflg)==28,'bwpt flag wrong');
        assert(length(bwftflg)==27,'bwft flag wrong');
        assert(length(bwcflg)==25, 'bwc flag wrong');
        assert(length(bgwptflg)==10,'bgwpt flag wrong');
        assert(length(bgwftflg)==10,'bgwft flag wrong');
        assert(length(bgwcflg)==10, 'bgwc flag wrong');
    end

    % in the below, there is a zero after the "sdemog" block in each alternative. That's because consumption gets added on separately at the end of the utility expression
    alpha     = b(1);
    galpha    = b(2);
    b2        = b(b2flg);
    b2temp    = [b2(1:sdemog+3);0;b2(sdemog+4:end-7);0;b2(end-6:end)];% first 0: consumption; second 0: white_collar
    b2tflg    = find(b2temp~=0);
    b4s       = b(b4sflg);
    b4stemp   = [b4s(1:sdemog+3);0;b4s(sdemog+4:end-7);0;b4s(end-6:end)];
    b4stflg   = find(b4stemp~=0);
    b4ns      = b(b4nsflg);
    b4nstemp  = [b4ns(1:sdemog+3);0;b4ns(sdemog+4:end-7);0;b4ns(end-6:end)];
    b4nstflg  = find(b4nstemp~=0);
    bwpt      = b(bwptflg);
    bwpttemp  = [bwpt(1:sdemog);0;bwpt(sdemog+1:sdemog+2);0;bwpt(sdemog+3:end-3);0;0;0;0;bwpt(end-2:end)];% first two 0s: academic ability, consumption; second four 0s: workPT, workFT, workPT*white_collar, workFT*white_collar
    bwpttflg  = find(bwpttemp~=0);
    bwft      = b(bwftflg);
    bwfttemp  = [bwft(1:sdemog);0;bwft(sdemog+1:sdemog+2);0;bwft(sdemog+3:end-3);0;0;0;0;0;bwft(end-2:end)];% first two 0s: academic ability, consumption; second 5 0s: white_collar, workPT, workFT, wPT*WC, wFT*WC
    bwfttflg  = find(bwfttemp~=0);
    bwc       = b(bwcflg);
    bwctemp   = [bwc(1:sdemog);0;0;0;0;bwc(sdemog+1:end-3);0;0;0;0;0;bwc(end-2:end)];% first four 0s: academic ability, accum. debt (and squared), consumption; second 5 0s: white_collar, workPT, workFT, wPT*WC, wFT*WC 
    bwctflg   = find(bwctemp~=0);

    bgwpt     = b(bgwptflg);
    bgwpttemp = bgwpt(1:10);
    bgwpttflg = find(bgwpttemp~=0);
    bgwft     = b(bgwftflg);
    bgwfttemp = bgwft(1:10);
    bgwfttflg = find(bgwfttemp~=0);
    bgwc      = b(bgwcflg);
    bgwctemp  = bgwc(1:10);
    bgwctflg  = find(bgwctemp~=0);

    %disp('bgwfttflg');
    %disp(bgwfttflg);
    %assert(all(bgwfttflg(:)==[1:10]'))
    %assert(all(bgwpttflg(:)==[1:10]'))
    %assert(all(bgwctflg(:)==[1:10]'))

    % indices to flag consumption and non-consumption covariates in matrices
    cidx           = sdemog+4;
    gcidx          = 11;

    % "utility" for 2-year college:
    vng2ftbc = (ut.X2ftbc(flag,:)*(b2temp+bwfttemp        )+ut.X2ftbc(flag,cidx)*alpha);
    vng2ftwc = (ut.X2ftwc(flag,:)*(b2temp+bwfttemp+bwctemp)+ut.X2ftwc(flag,cidx)*alpha);
    vng2ptbc = (ut.X2ptbc(flag,:)*(b2temp+bwpttemp        )+ut.X2ptbc(flag,cidx)*alpha);
    vng2ptwc = (ut.X2ptwc(flag,:)*(b2temp+bwpttemp+bwctemp)+ut.X2ptwc(flag,cidx)*alpha);
    vng2     = (ut.X2nw(  flag,:)*(b2temp                 )+ut.X2nw(  flag,cidx)*alpha);

    % "utility" for 4-year college: Science majors
    vng4sftbc = (ut.X4sftbc(flag,:)*(b4stemp+bwfttemp        )+ut.X4sftbc(flag,cidx)*alpha);
    vng4sftwc = (ut.X4sftwc(flag,:)*(b4stemp+bwfttemp+bwctemp)+ut.X4sftwc(flag,cidx)*alpha);
    vng4sptbc = (ut.X4sptbc(flag,:)*(b4stemp+bwpttemp        )+ut.X4sptbc(flag,cidx)*alpha);
    vng4sptwc = (ut.X4sptwc(flag,:)*(b4stemp+bwpttemp+bwctemp)+ut.X4sptwc(flag,cidx)*alpha);
    vng4s     = (ut.X4snw(  flag,:)*(b4stemp                 )+ut.X4snw(  flag,cidx)*alpha);

    % Non-Science majors
    vng4nsftbc = (ut.X4nsftbc(flag,:)*(b4nstemp+bwfttemp        )+ut.X4nsftbc(flag,cidx)*alpha);
    vng4nsftwc = (ut.X4nsftwc(flag,:)*(b4nstemp+bwfttemp+bwctemp)+ut.X4nsftwc(flag,cidx)*alpha);
    vng4nsptbc = (ut.X4nsptbc(flag,:)*(b4nstemp+bwpttemp        )+ut.X4nsptbc(flag,cidx)*alpha);
    vng4nsptwc = (ut.X4nsptwc(flag,:)*(b4nstemp+bwpttemp+bwctemp)+ut.X4nsptwc(flag,cidx)*alpha);
    vng4ns     = (ut.X4nsnw(  flag,:)*(b4nstemp                 )+ut.X4nsnw(  flag,cidx)*alpha);

    % "utility" for working
    vngwptbc = (ut.Xngwptbc(flag,:)*(bwpttemp        )+ut.Xngwptbc(flag,cidx)*alpha);
    vngwptwc = (ut.Xngwptwc(flag,:)*(bwpttemp+bwctemp)+ut.Xngwptwc(flag,cidx)*alpha);
    vngwftbc = (ut.Xngwftbc(flag,:)*(bwfttemp        )+ut.Xngwftbc(flag,cidx)*alpha);
    vngwftwc = (ut.Xngwftwc(flag,:)*(bwfttemp+bwctemp)+ut.Xngwftwc(flag,cidx)*alpha);

    % "utility" for graduate options
    vgwptbc = (ut.Xngwptbc( flag,:)*(bwpttemp        )+ut.Xgwptbc( flag,bgwpttflg)*(bgwpttemp         )+ut.Xgwptbc( flag,gcidx)*galpha);
    vgwptwc = (ut.Xngwptwc( flag,:)*(bwpttemp+bwctemp)+ut.Xgwptwc( flag,bgwpttflg)*(bgwpttemp+bgwctemp)+ut.Xgwptwc( flag,gcidx)*galpha);
    vgwftbc = (ut.Xngwftbc( flag,:)*(bwfttemp        )+ut.Xgwftbc( flag,bgwfttflg)*(bgwfttemp         )+ut.Xgwftbc( flag,gcidx)*galpha);
    vgwftwc = (ut.Xngwftwc( flag,:)*(bwfttemp+bwctemp)+ut.Xgwftwc( flag,bgwfttflg)*(bgwfttemp+bgwctemp)+ut.Xgwftwc( flag,gcidx)*galpha);

    %%% debugging:
    %latidx = Xo(flag & prev_WC==0,:)*bo;
    %allmatNG = [vng2ftbc vng2ftwc vng2ptbc vng2ptwc vng2 vng4sftbc vng4sftwc vng4sptbc vng4sptwc vng4s vng4nsftbc vng4nsftwc vng4nsptbc vng4nsptwc vng4ns vngwptbc vngwptwc vngwftbc vngwftwc];
    %allmatG = [vgwptbc vgwptwc vgwftbc vgwftwc];
    %disp('sum stats non-grad util');
    %sumopt = struct('Weights',q(flag).*(grad_4yr(flag)==0));
    %summarize(allmatNG,sumopt);
    %disp('sum stats grad util');
    %sumopt = struct('Weights',q(flag).*(grad_4yr(flag)==1));
    %summarize(allmatG,sumopt);
    %disp('sum stats all util');
    %sumopt = struct('Weights',q(flag));
    %summarize([allmatNG allmatG],sumopt);
    %disp('sum stats lambda index');
    %sumopt = struct('Weights',q(flag & prev_WC==0));
    %summarize(latidx,sumopt);
    %disp('offer parameters');
    %bo
    %disp('choice parameters');
    %b
    %blablabla

    %% Form likelihood (3 parts)
    % Part 1a
    log_dem_ngno = log(1+exp(vng2ftbc)+exp(vng2ptbc)+exp(vng2)+exp(vng4sftbc)+exp(vng4sptbc)+exp(vng4s)+exp(vng4nsftbc)+exp(vng4nsptbc)+exp(vng4ns)+exp(vngwptbc)+exp(vngwftbc));
    log_p2ftbc_ngno   = vng2ftbc  -log_dem_ngno;
    log_p2ptbc_ngno   = vng2ptbc  -log_dem_ngno;
    log_p2_ngno       = vng2      -log_dem_ngno;
    log_p4sftbc_ngno  = vng4sftbc -log_dem_ngno;
    log_p4sptbc_ngno  = vng4sptbc -log_dem_ngno;
    log_p4s_ngno      = vng4s     -log_dem_ngno;
    log_p4nsftbc_ngno = vng4nsftbc-log_dem_ngno;
    log_p4nsptbc_ngno = vng4nsptbc-log_dem_ngno;
    log_p4ns_ngno     = vng4ns    -log_dem_ngno;
    log_pwptbc_ngno   = vngwptbc  -log_dem_ngno;
    log_pwftbc_ngno   = vngwftbc  -log_dem_ngno;
    log_ph_ngno       =           -log_dem_ngno;

    % Part 1b
    log_dem_ngof = log(1+exp(vng2ftbc)+exp(vng2ftwc)+exp(vng2ptbc)+exp(vng2ptwc)+exp(vng2)+exp(vng4sftbc)+exp(vng4sftwc)+exp(vng4sptbc)+exp(vng4sptwc)+exp(vng4s)+exp(vng4nsftbc)+exp(vng4nsftwc)+exp(vng4nsptbc)+exp(vng4nsptwc)+exp(vng4ns)+exp(vngwptbc)+exp(vngwptwc)+exp(vngwftbc)+exp(vngwftwc));
    log_p2ftbc_ngof   = vng2ftbc  -log_dem_ngof;
    log_p2ftwc_ngof   = vng2ftwc  -log_dem_ngof;
    log_p2ptbc_ngof   = vng2ptbc  -log_dem_ngof;
    log_p2ptwc_ngof   = vng2ptwc  -log_dem_ngof;
    log_p2_ngof       = vng2      -log_dem_ngof;
    log_p4sftbc_ngof  = vng4sftbc -log_dem_ngof;
    log_p4sftwc_ngof  = vng4sftwc -log_dem_ngof;
    log_p4sptbc_ngof  = vng4sptbc -log_dem_ngof;
    log_p4sptwc_ngof  = vng4sptwc -log_dem_ngof;
    log_p4s_ngof      = vng4s     -log_dem_ngof;
    log_p4nsftbc_ngof = vng4nsftbc-log_dem_ngof;
    log_p4nsftwc_ngof = vng4nsftwc-log_dem_ngof;
    log_p4nsptbc_ngof = vng4nsptbc-log_dem_ngof;
    log_p4nsptwc_ngof = vng4nsptwc-log_dem_ngof;
    log_p4ns_ngof     = vng4ns    -log_dem_ngof;
    log_pwptbc_ngof   = vngwptbc  -log_dem_ngof;
    log_pwptwc_ngof   = vngwptwc  -log_dem_ngof;
    log_pwftbc_ngof   = vngwftbc  -log_dem_ngof;
    log_pwftwc_ngof   = vngwftwc  -log_dem_ngof;
    log_ph_ngof       =           -log_dem_ngof;

    % Part 2a
    log_dem_grno     = log(1+exp(vgwptbc)+exp(vgwftbc));
    log_pwptbc_grno  = vgwptbc-log_dem_grno;
    log_pwftbc_grno  = vgwftbc-log_dem_grno;
    log_ph_grno      =        -log_dem_grno;

    % Part 2b
    log_dem_grof     = log(1+exp(vgwptbc)+exp(vgwptwc)+exp(vgwftbc)+exp(vgwftwc));
    log_pwptbc_grof  = vgwptbc-log_dem_grof;
    log_pwptwc_grof  = vgwptwc-log_dem_grof;
    log_pwftbc_grof  = vgwftbc-log_dem_grof;
    log_pwftwc_grof  = vgwftwc-log_dem_grof;
    log_ph_grof      =        -log_dem_grof;

    log_like_ngno = (Y(flag)==1).*log_p2ftbc_ngno    +...
                    (Y(flag)==3).*log_p2ptbc_ngno    +...
                    (Y(flag)==5).*log_p2_ngno        +...
                    (Y(flag)==6).*log_p4sftbc_ngno   +...
                    (Y(flag)==8).*log_p4sptbc_ngno   +...
                    (Y(flag)==10).*log_p4s_ngno      +...
                    (Y(flag)==11).*log_p4nsftbc_ngno +...
                    (Y(flag)==13).*log_p4nsptbc_ngno +...
                    (Y(flag)==15).*log_p4ns_ngno     +...
                    (Y(flag)==16).*log_pwptbc_ngno   +...
                    (Y(flag)==18).*log_pwftbc_ngno   +...
                    (Y(flag)==20).*log_ph_ngno;

    log_like_ngof = (Y(flag)==1).*log_p2ftbc_ngof    +...
                    (Y(flag)==2).*log_p2ftwc_ngof    +...
                    (Y(flag)==3).*log_p2ptbc_ngof    +...
                    (Y(flag)==4).*log_p2ptwc_ngof    +...
                    (Y(flag)==5).*log_p2_ngof        +...
                    (Y(flag)==6).*log_p4sftbc_ngof   +...
                    (Y(flag)==7).*log_p4sftwc_ngof   +...
                    (Y(flag)==8).*log_p4sptbc_ngof   +...
                    (Y(flag)==9).*log_p4sptwc_ngof   +...
                    (Y(flag)==10).*log_p4s_ngof      +...
                    (Y(flag)==11).*log_p4nsftbc_ngof +...
                    (Y(flag)==12).*log_p4nsftwc_ngof +...
                    (Y(flag)==13).*log_p4nsptbc_ngof +...
                    (Y(flag)==14).*log_p4nsptwc_ngof +...
                    (Y(flag)==15).*log_p4ns_ngof     +...
                    (Y(flag)==16).*log_pwptbc_ngof   +...
                    (Y(flag)==17).*log_pwptwc_ngof   +...
                    (Y(flag)==18).*log_pwftbc_ngof   +...
                    (Y(flag)==19).*log_pwftwc_ngof   +...
                    (Y(flag)==20).*log_ph_ngof;

    log_like_grno = (Y(flag)==16).*log_pwptbc_grno  +...
                    (Y(flag)==18).*log_pwftbc_grno  +...
                    (Y(flag)==20).*log_ph_grno;

    log_like_grof = (Y(flag)==16).*log_pwptbc_grof  +...
                    (Y(flag)==17).*log_pwptwc_grof  +...
                    (Y(flag)==18).*log_pwftbc_grof  +...
                    (Y(flag)==19).*log_pwftwc_grof  +...
                    (Y(flag)==20).*log_ph_grof;

    like_o  = exp(log_like_ngof.*(grad_4yr(flag)==0)+log_like_grof.*(grad_4yr(flag)==1)); % equation G.4 in appendix
    like_no = exp(log_like_ngno.*(grad_4yr(flag)==0)+log_like_grno.*(grad_4yr(flag)==1)); % equation G.4 in appendix

    ll0     = log(lambda(flag).*like_o + (1-lambda(flag)).*like_no);
    ll1     = log(lambda(flag).*like_o);

    like = -ctranspose(q(flag))*( ( ll0.*(in_white_collar(flag)==0 & prev_WC(flag)==0) ) + ( ll1.*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) );

    % Analytical gradient
    grad = zeros(size(bjoint));

    dem_no  = lambda(flag).*like_o + (1-lambda(flag)).*like_no;

    % create matrix of terms that hold the q*(d-P)*grad terms that always show up in the gradient
    temp = zeros(sum(flag),20,4);
    temp(:,1 ,1) = ( q(flag).*((Y(flag)==1 )-exp(log_p2ftbc_ngno  )).*(grad_4yr(flag)==0) );
    temp(:,3 ,1) = ( q(flag).*((Y(flag)==3 )-exp(log_p2ptbc_ngno  )).*(grad_4yr(flag)==0) );
    temp(:,5 ,1) = ( q(flag).*((Y(flag)==5 )-exp(log_p2_ngno      )).*(grad_4yr(flag)==0) );
    temp(:,6 ,1) = ( q(flag).*((Y(flag)==6 )-exp(log_p4sftbc_ngno )).*(grad_4yr(flag)==0) );
    temp(:,8 ,1) = ( q(flag).*((Y(flag)==8 )-exp(log_p4sptbc_ngno )).*(grad_4yr(flag)==0) );
    temp(:,10,1) = ( q(flag).*((Y(flag)==10)-exp(log_p4s_ngno     )).*(grad_4yr(flag)==0) );
    temp(:,11,1) = ( q(flag).*((Y(flag)==11)-exp(log_p4nsftbc_ngno)).*(grad_4yr(flag)==0) );
    temp(:,13,1) = ( q(flag).*((Y(flag)==13)-exp(log_p4nsptbc_ngno)).*(grad_4yr(flag)==0) );
    temp(:,15,1) = ( q(flag).*((Y(flag)==15)-exp(log_p4ns_ngno    )).*(grad_4yr(flag)==0) );
    temp(:,16,1) = ( q(flag).*((Y(flag)==16)-exp(log_pwptbc_ngno  )).*(grad_4yr(flag)==0) );
    temp(:,18,1) = ( q(flag).*((Y(flag)==18)-exp(log_pwftbc_ngno  )).*(grad_4yr(flag)==0) );
    temp(:,20,1) = ( q(flag).*((Y(flag)==20)-exp(log_ph_ngno      )).*(grad_4yr(flag)==0) );

    temp(:,1 ,2) = ( q(flag).*((Y(flag)==1 )-exp(log_p2ftbc_ngof  )).*(grad_4yr(flag)==0) );
    temp(:,2 ,2) = ( q(flag).*((Y(flag)==2 )-exp(log_p2ftwc_ngof  )).*(grad_4yr(flag)==0) );
    temp(:,3 ,2) = ( q(flag).*((Y(flag)==3 )-exp(log_p2ptbc_ngof  )).*(grad_4yr(flag)==0) );
    temp(:,4 ,2) = ( q(flag).*((Y(flag)==4 )-exp(log_p2ptwc_ngof  )).*(grad_4yr(flag)==0) );
    temp(:,5 ,2) = ( q(flag).*((Y(flag)==5 )-exp(log_p2_ngof      )).*(grad_4yr(flag)==0) );
    temp(:,6 ,2) = ( q(flag).*((Y(flag)==6 )-exp(log_p4sftbc_ngof )).*(grad_4yr(flag)==0) );
    temp(:,7 ,2) = ( q(flag).*((Y(flag)==7 )-exp(log_p4sftwc_ngof )).*(grad_4yr(flag)==0) );
    temp(:,8 ,2) = ( q(flag).*((Y(flag)==8 )-exp(log_p4sptbc_ngof )).*(grad_4yr(flag)==0) );
    temp(:,9 ,2) = ( q(flag).*((Y(flag)==9 )-exp(log_p4sptwc_ngof )).*(grad_4yr(flag)==0) );
    temp(:,10,2) = ( q(flag).*((Y(flag)==10)-exp(log_p4s_ngof     )).*(grad_4yr(flag)==0) );
    temp(:,11,2) = ( q(flag).*((Y(flag)==11)-exp(log_p4nsftbc_ngof)).*(grad_4yr(flag)==0) );
    temp(:,12,2) = ( q(flag).*((Y(flag)==12)-exp(log_p4nsftwc_ngof)).*(grad_4yr(flag)==0) );
    temp(:,13,2) = ( q(flag).*((Y(flag)==13)-exp(log_p4nsptbc_ngof)).*(grad_4yr(flag)==0) );
    temp(:,14,2) = ( q(flag).*((Y(flag)==14)-exp(log_p4nsptwc_ngof)).*(grad_4yr(flag)==0) );
    temp(:,15,2) = ( q(flag).*((Y(flag)==15)-exp(log_p4ns_ngof    )).*(grad_4yr(flag)==0) );
    temp(:,16,2) = ( q(flag).*((Y(flag)==16)-exp(log_pwptbc_ngof  )).*(grad_4yr(flag)==0) );
    temp(:,17,2) = ( q(flag).*((Y(flag)==17)-exp(log_pwptwc_ngof  )).*(grad_4yr(flag)==0) );
    temp(:,18,2) = ( q(flag).*((Y(flag)==18)-exp(log_pwftbc_ngof  )).*(grad_4yr(flag)==0) );
    temp(:,19,2) = ( q(flag).*((Y(flag)==19)-exp(log_pwftwc_ngof  )).*(grad_4yr(flag)==0) );
    temp(:,20,2) = ( q(flag).*((Y(flag)==20)-exp(log_ph_ngof      )).*(grad_4yr(flag)==0) );

    temp(:,16,3) = ( q(flag).*((Y(flag)==16)-exp(log_pwptbc_grno  )).*(grad_4yr(flag)==1) );
    temp(:,18,3) = ( q(flag).*((Y(flag)==18)-exp(log_pwftbc_grno  )).*(grad_4yr(flag)==1) );
    temp(:,20,3) = ( q(flag).*((Y(flag)==20)-exp(log_ph_grno      )).*(grad_4yr(flag)==1) );

    temp(:,16,4) = ( q(flag).*((Y(flag)==16)-exp(log_pwptbc_grof  )).*(grad_4yr(flag)==1) );
    temp(:,17,4) = ( q(flag).*((Y(flag)==17)-exp(log_pwptwc_grof  )).*(grad_4yr(flag)==1) );
    temp(:,18,4) = ( q(flag).*((Y(flag)==18)-exp(log_pwftbc_grof  )).*(grad_4yr(flag)==1) );
    temp(:,19,4) = ( q(flag).*((Y(flag)==19)-exp(log_pwftwc_grof  )).*(grad_4yr(flag)==1) );
    temp(:,20,4) = ( q(flag).*((Y(flag)==20)-exp(log_ph_grof      )).*(grad_4yr(flag)==1) );

    grad_ngno=zeros(size(b));
    grad_ngno(1,1)       = -(ut.X2ftbc  (flag,cidx    ))'*( temp(:, 1,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X2ptbc  (flag,cidx    ))'*( temp(:, 3,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X2nw    (flag,cidx    ))'*( temp(:, 5,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4sftbc (flag,cidx    ))'*( temp(:, 6,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4sptbc (flag,cidx    ))'*( temp(:, 8,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4snw   (flag,cidx    ))'*( temp(:,10,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4nsftbc(flag,cidx    ))'*( temp(:,11,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4nsptbc(flag,cidx    ))'*( temp(:,13,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4nsnw  (flag,cidx    ))'*( temp(:,15,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.Xngwptbc(flag,cidx    ))'*( temp(:,16,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.Xngwftbc(flag,cidx    ))'*( temp(:,18,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no);
    grad_ngno(b2flg,1)   = -(ut.X2ftbc  (flag,b2tflg  ))'*( temp(:, 1,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X2ptbc  (flag,b2tflg  ))'*( temp(:, 3,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X2nw    (flag,b2tflg  ))'*( temp(:, 5,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no);
    grad_ngno(b4sflg,1)  = -(ut.X4sftbc (flag,b4stflg ))'*( temp(:, 6,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4sptbc (flag,b4stflg ))'*( temp(:, 8,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4snw   (flag,b4stflg ))'*( temp(:,10,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no);
    grad_ngno(b4nsflg,1) = -(ut.X4nsftbc(flag,b4nstflg))'*( temp(:,11,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4nsptbc(flag,b4nstflg))'*( temp(:,13,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4nsnw  (flag,b4nstflg))'*( temp(:,15,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no);
    grad_ngno(bwptflg,1) = -(ut.X2ptbc  (flag,bwpttflg))'*( temp(:, 3,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4sptbc (flag,bwpttflg))'*( temp(:, 8,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4nsptbc(flag,bwpttflg))'*( temp(:,13,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.Xngwptbc(flag,bwpttflg))'*( temp(:,16,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no);
    grad_ngno(bwftflg,1) = -(ut.X2ftbc  (flag,bwfttflg))'*( temp(:, 1,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4sftbc (flag,bwfttflg))'*( temp(:, 6,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.X4nsftbc(flag,bwfttflg))'*( temp(:,11,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no) - ...
                            (ut.Xngwftbc(flag,bwfttflg))'*( temp(:,18,1).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no);

    grad_ngof_o=zeros(size(b));
    grad_ngof_o(1,1)       = -(ut.X2ftbc  (flag,cidx    ))'*( temp(:, 1,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ftwc  (flag,cidx    ))'*( temp(:, 2,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ptbc  (flag,cidx    ))'*( temp(:, 3,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ptwc  (flag,cidx    ))'*( temp(:, 4,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2nw    (flag,cidx    ))'*( temp(:, 5,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sftbc (flag,cidx    ))'*( temp(:, 6,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sftwc (flag,cidx    ))'*( temp(:, 7,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sptbc (flag,cidx    ))'*( temp(:, 8,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sptwc (flag,cidx    ))'*( temp(:, 9,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4snw   (flag,cidx    ))'*( temp(:,10,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsftbc(flag,cidx    ))'*( temp(:,11,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsftwc(flag,cidx    ))'*( temp(:,12,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsptbc(flag,cidx    ))'*( temp(:,13,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsptwc(flag,cidx    ))'*( temp(:,14,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsnw  (flag,cidx    ))'*( temp(:,15,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwptbc(flag,cidx    ))'*( temp(:,16,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwptwc(flag,cidx    ))'*( temp(:,17,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwftbc(flag,cidx    ))'*( temp(:,18,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwftwc(flag,cidx    ))'*( temp(:,19,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_ngof_o(b2flg,1)   = -(ut.X2ftbc  (flag,b2tflg  ))'*( temp(:, 1,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ftwc  (flag,b2tflg  ))'*( temp(:, 2,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ptbc  (flag,b2tflg  ))'*( temp(:, 3,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ptwc  (flag,b2tflg  ))'*( temp(:, 4,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2nw    (flag,b2tflg  ))'*( temp(:, 5,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_ngof_o(b4sflg,1)  = -(ut.X4sftbc (flag,b4stflg ))'*( temp(:, 6,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sftwc (flag,b4stflg ))'*( temp(:, 7,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sptbc (flag,b4stflg ))'*( temp(:, 8,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sptwc (flag,b4stflg ))'*( temp(:, 9,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4snw   (flag,b4stflg ))'*( temp(:,10,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_ngof_o(b4nsflg,1) = -(ut.X4nsftbc(flag,b4nstflg))'*( temp(:,11,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsftwc(flag,b4nstflg))'*( temp(:,12,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsptbc(flag,b4nstflg))'*( temp(:,13,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsptwc(flag,b4nstflg))'*( temp(:,14,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsnw  (flag,b4nstflg))'*( temp(:,15,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_ngof_o(bwptflg,1) = -(ut.X2ptbc  (flag,bwpttflg))'*( temp(:, 3,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ptwc  (flag,bwpttflg))'*( temp(:, 4,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sptbc (flag,bwpttflg))'*( temp(:, 8,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sptwc (flag,bwpttflg))'*( temp(:, 9,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsptbc(flag,bwpttflg))'*( temp(:,13,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsptwc(flag,bwpttflg))'*( temp(:,14,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwptbc(flag,bwpttflg))'*( temp(:,16,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwptwc(flag,bwpttflg))'*( temp(:,17,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_ngof_o(bwftflg,1) = -(ut.X2ftbc  (flag,bwfttflg))'*( temp(:, 1,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ftwc  (flag,bwfttflg))'*( temp(:, 2,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sftbc (flag,bwfttflg))'*( temp(:, 6,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sftwc (flag,bwfttflg))'*( temp(:, 7,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsftbc(flag,bwfttflg))'*( temp(:,11,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsftwc(flag,bwfttflg))'*( temp(:,12,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwftbc(flag,bwfttflg))'*( temp(:,18,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwftwc(flag,bwfttflg))'*( temp(:,19,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_ngof_o(bwcflg,1)  = -(ut.X2ftwc  (flag,bwctflg ))'*( temp(:, 2,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X2ptwc  (flag,bwctflg ))'*( temp(:, 4,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sftwc (flag,bwctflg ))'*( temp(:, 7,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4sptwc (flag,bwctflg ))'*( temp(:, 9,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsftwc(flag,bwctflg ))'*( temp(:,12,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.X4nsptwc(flag,bwctflg ))'*( temp(:,14,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwptwc(flag,bwctflg ))'*( temp(:,17,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                              (ut.Xngwftwc(flag,bwctflg ))'*( temp(:,19,2).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );

    grad_ngof_no=zeros(size(b));
    grad_ngof_no(1,1)       = -(ut.X2ftbc  (flag,cidx    ))'*( temp(:, 1,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ftwc  (flag,cidx    ))'*( temp(:, 2,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ptbc  (flag,cidx    ))'*( temp(:, 3,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ptwc  (flag,cidx    ))'*( temp(:, 4,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2nw    (flag,cidx    ))'*( temp(:, 5,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sftbc (flag,cidx    ))'*( temp(:, 6,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sftwc (flag,cidx    ))'*( temp(:, 7,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sptbc (flag,cidx    ))'*( temp(:, 8,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sptwc (flag,cidx    ))'*( temp(:, 9,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4snw   (flag,cidx    ))'*( temp(:,10,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsftbc(flag,cidx    ))'*( temp(:,11,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsftwc(flag,cidx    ))'*( temp(:,12,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsptbc(flag,cidx    ))'*( temp(:,13,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsptwc(flag,cidx    ))'*( temp(:,14,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsnw  (flag,cidx    ))'*( temp(:,15,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwptbc(flag,cidx    ))'*( temp(:,16,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwptwc(flag,cidx    ))'*( temp(:,17,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwftbc(flag,cidx    ))'*( temp(:,18,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwftwc(flag,cidx    ))'*( temp(:,19,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_ngof_no(b2flg,1)   = -(ut.X2ftbc  (flag,b2tflg  ))'*( temp(:, 1,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ftwc  (flag,b2tflg  ))'*( temp(:, 2,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ptbc  (flag,b2tflg  ))'*( temp(:, 3,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ptwc  (flag,b2tflg  ))'*( temp(:, 4,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2nw    (flag,b2tflg  ))'*( temp(:, 5,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_ngof_no(b4sflg,1)  = -(ut.X4sftbc (flag,b4stflg ))'*( temp(:, 6,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sftwc (flag,b4stflg ))'*( temp(:, 7,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sptbc (flag,b4stflg ))'*( temp(:, 8,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sptwc (flag,b4stflg ))'*( temp(:, 9,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4snw   (flag,b4stflg ))'*( temp(:,10,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_ngof_no(b4nsflg,1) = -(ut.X4nsftbc(flag,b4nstflg))'*( temp(:,11,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsftwc(flag,b4nstflg))'*( temp(:,12,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsptbc(flag,b4nstflg))'*( temp(:,13,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsptwc(flag,b4nstflg))'*( temp(:,14,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsnw  (flag,b4nstflg))'*( temp(:,15,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_ngof_no(bwptflg,1) = -(ut.X2ptbc  (flag,bwpttflg))'*( temp(:, 3,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ptwc  (flag,bwpttflg))'*( temp(:, 4,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sptbc (flag,bwpttflg))'*( temp(:, 8,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sptwc (flag,bwpttflg))'*( temp(:, 9,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsptbc(flag,bwpttflg))'*( temp(:,13,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsptwc(flag,bwpttflg))'*( temp(:,14,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwptbc(flag,bwpttflg))'*( temp(:,16,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwptwc(flag,bwpttflg))'*( temp(:,17,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_ngof_no(bwftflg,1) = -(ut.X2ftbc  (flag,bwfttflg))'*( temp(:, 1,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ftwc  (flag,bwfttflg))'*( temp(:, 2,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sftbc (flag,bwfttflg))'*( temp(:, 6,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sftwc (flag,bwfttflg))'*( temp(:, 7,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsftbc(flag,bwfttflg))'*( temp(:,11,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsftwc(flag,bwfttflg))'*( temp(:,12,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwftbc(flag,bwfttflg))'*( temp(:,18,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwftwc(flag,bwfttflg))'*( temp(:,19,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_ngof_no(bwcflg,1)  = -(ut.X2ftwc  (flag,bwctflg ))'*( temp(:, 2,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X2ptwc  (flag,bwctflg ))'*( temp(:, 4,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sftwc (flag,bwctflg ))'*( temp(:, 7,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4sptwc (flag,bwctflg ))'*( temp(:, 9,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsftwc(flag,bwctflg ))'*( temp(:,12,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.X4nsptwc(flag,bwctflg ))'*( temp(:,14,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwptwc(flag,bwctflg ))'*( temp(:,17,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                               (ut.Xngwftwc(flag,bwctflg ))'*( temp(:,19,2).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );

    grad_grno=zeros(size(b));
    grad_grno(2,1)        = -(ut.Xgwptbc (flag,gcidx    ))'*( temp(:,16,3).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no ) - ...
                             (ut.Xgwftbc (flag,gcidx    ))'*( temp(:,18,3).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no );
    grad_grno(bwptflg,1)  = -(ut.Xngwptbc(flag,bwpttflg ))'*( temp(:,16,3).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no );
    grad_grno(bwftflg,1)  = -(ut.Xngwftbc(flag,bwfttflg ))'*( temp(:,18,3).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no );
    grad_grno(bgwptflg,1) = -(ut.Xgwptbc (flag,bgwpttflg))'*( temp(:,16,3).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no );
    grad_grno(bgwftflg,1) = -(ut.Xgwftbc (flag,bgwfttflg))'*( temp(:,18,3).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*(1-lambda(flag)).*like_no./dem_no );

    grad_grof_no=zeros(size(b));
    grad_grof_no(2,1)        = -(ut.Xgwptbc (flag,gcidx    ))'*( temp(:,16,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                                (ut.Xgwptwc (flag,gcidx    ))'*( temp(:,17,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                                (ut.Xgwftbc (flag,gcidx    ))'*( temp(:,18,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                                (ut.Xgwftwc (flag,gcidx    ))'*( temp(:,19,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_grof_no(bwptflg,1)  = -(ut.Xngwptbc(flag,bwpttflg ))'*( temp(:,16,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                                (ut.Xngwptwc(flag,bwpttflg ))'*( temp(:,17,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_grof_no(bwftflg,1)  = -(ut.Xngwftbc(flag,bwfttflg ))'*( temp(:,18,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                                (ut.Xngwftwc(flag,bwfttflg ))'*( temp(:,19,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_grof_no(bwcflg,1)   = -(ut.Xngwptwc(flag,bwctflg  ))'*( temp(:,17,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                                (ut.Xngwftwc(flag,bwctflg  ))'*( temp(:,19,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_grof_no(bgwptflg,1) = -(ut.Xgwptbc (flag,bgwpttflg))'*( temp(:,16,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                                (ut.Xgwptwc (flag,bgwpttflg))'*( temp(:,17,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_grof_no(bgwftflg,1) = -(ut.Xgwftbc (flag,bgwfttflg))'*( temp(:,18,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                                (ut.Xgwftwc (flag,bgwfttflg))'*( temp(:,19,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );
    grad_grof_no(bgwcflg,1)  = -(ut.Xgwptwc (flag,bgwctflg ))'*( temp(:,17,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no ) - ...
                                (ut.Xgwftwc (flag,bgwctflg ))'*( temp(:,19,4).*(in_white_collar(flag)==0 & prev_WC(flag)==0).*lambda(flag).*like_o./dem_no );

    grad_grof_o=zeros(size(b));
    grad_grof_o(2,1)        = -(ut.Xgwptbc (flag,gcidx    ))'*( temp(:,16,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                               (ut.Xgwptwc (flag,gcidx    ))'*( temp(:,17,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                               (ut.Xgwftbc (flag,gcidx    ))'*( temp(:,18,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                               (ut.Xgwftwc (flag,gcidx    ))'*( temp(:,19,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_grof_o(bwptflg,1)  = -(ut.Xngwptbc(flag,bwpttflg ))'*( temp(:,16,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                               (ut.Xngwptwc(flag,bwpttflg ))'*( temp(:,17,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_grof_o(bwftflg,1)  = -(ut.Xngwftbc(flag,bwfttflg ))'*( temp(:,18,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                               (ut.Xngwftwc(flag,bwfttflg ))'*( temp(:,19,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_grof_o(bwcflg,1)   = -(ut.Xngwptwc(flag,bwctflg  ))'*( temp(:,17,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                               (ut.Xngwftwc(flag,bwctflg  ))'*( temp(:,19,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_grof_o(bgwptflg,1) = -(ut.Xgwptbc (flag,bgwpttflg))'*( temp(:,16,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                               (ut.Xgwptwc (flag,bgwpttflg))'*( temp(:,17,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_grof_o(bgwftflg,1) = -(ut.Xgwftbc (flag,bgwfttflg))'*( temp(:,18,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                               (ut.Xgwftwc (flag,bgwfttflg))'*( temp(:,19,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );
    grad_grof_o(bgwcflg,1)  = -(ut.Xgwptwc (flag,bgwctflg ))'*( temp(:,17,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) ) - ...
                               (ut.Xgwftwc (flag,bgwctflg ))'*( temp(:,19,4).*(in_white_collar(flag)==1 | prev_WC(flag)==1) );

    grad_no = grad_ngno+grad_grno+grad_ngof_no+grad_grof_no;
    grad_o  = grad_ngof_o+grad_grof_o;
    grad(idx1) = grad_no+grad_o;
    grad(idx2) = -(Xo(flag & prev_WC==0,:)'*( q(flag & prev_WC==0).*(1-lambda(flag & prev_WC==0)).*(in_white_collar(flag & prev_WC==0)==1 | prev_WC(flag & prev_WC==0)==1) )) - (Xo(flag & prev_WC==0,:)'*( q(flag & prev_WC==0).*lambda(flag & prev_WC==0).*( (like_o(prev_WC(flag)==0)-like_no(prev_WC(flag)==0))./( like_o(prev_WC(flag)==0).*exp(Xo(flag & prev_WC==0,:)*bo)+like_no(prev_WC(flag)==0) ) ).*(in_white_collar(flag & prev_WC==0)==0 & prev_WC(flag & prev_WC==0)==0) ));
end  
